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Abstract 



The three-dimensional secular behavior of a system composed of a central star and two 
massive planets is modeled semi-analytically in the frame of the general three-body problem. 
The main dynamical features of the system are presented in geometrical pictures allowing us 
to investigate a large domain of the phase space of this problem without time-expensive nu- 
merical integrations of the equations of motion and without any restriction on the magnitude 
of the planetary eccentricities, inclinations and mutual distance. Several regimes of motion of 
the system are observed. With respect to the secular angle Azu, possible motions are circula- 
tions, oscillations (around and 180°), and high eccentricity/inclination librations in secular 
resonances. With respect to the arguments of pericenter, w\ and cj 2 , possible motions are direct 
circulation and high-inclination libration around ±90° in the Lidov-Kozai resonance. The re- 
gions of transition between domains of different regimes of motion are characterized by chaotic 
behavior. 

We apply the analysis to the case of the two outer planets of the v Andromedae system, 
observed edge-on. The topology of the 3-D phase space of this system is investigated in detail 
by means of surfaces of section, periodic orbits and dynamical spectra, mapping techniques and 
numerical simulations. We obtain the general structure of the phase space, and the boundaries 
of the spatial secular stability. We find that this system is secularly stable in a large domain of 
eccentricities and inclinations. 



Key words: Celestial mechanics; Extrasolar planets; Planetary dynamics; Resonances 
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1 Introduction 



The discovery of extra-solar planets has brought into focus a long standing, but yet unsolved 
problem in Celestial Mechanics, namely the long-term stability of planetary systems. In fact, the 
dynamical evolution of the Solar System has been the object of exhaustive studies for at least two 
centuries, but its stability is not yet completely understood. The extra-solar systems provoke our 
imagination by their unusual configurations, which include Jupiter-mass planets at very small semi- 
major axes and unexpectedly large eccentricities, and raise new questions in our understanding of 
their dynamical stability. 

The first studies on the long-term stability of extra-solar planetary systems have been done by 
means of direct numerical integrations of the full equations of motion. In the last years, a large 
number of papers were published applying such approach to several planetary systems, in order 
to investigate their dynamical aspects. The numerical integration over very-long time intervals 
is a powerful approach for exploring the dynamical features of the systems with two and more 
planets. However, the correct interpretation of the results of numerical investigations is not simple 
and should be always founded on the basic concepts provided by analytical or semi-analytical 
investigations. Only these approaches are able to explore the whole parameter space of the problem 
and provide the decisive answer on the question of the global stability of the system under study. 

Until recently, the analytical models founded on the Lagrange-Laplace secular perturbation 
theory were commonly used in the study of the secular three-dimensional dynamics of the general 
three-body problem. These models are based on the classical expansion of the disturbing function 
in powers of eccentricities and inclinations, therefore their application is unsurprisingly restricted 
to the planetary systems with small eccentricities and inclinations (Murray and Dermott 1999). 
Even when a large number of the terms in the disturbing function is retained (Laskar 1987, Laskar 
and Robutel 1995), the convergence of the classical expansion remains limited to small values of the 
eccentricities, in accordance with the Sundman criterion (see Ferraz-Mello 1994), which establishes 
the convergence domains in the space of eccentricities for a given ratio of the semi- major axes. 

Analytical approaches based on the expansion of the disturbing function in a power series of 
the ratio of the planetary semi-major axes, a = a±/a2, have been used for triple star systems 
(Harrington 1969, Ford et al. 2000). The expansion to order a 3 , called octopole approximation, 
has been recently applied to study the secular dynamics of systems with two planets (Lee and Peale 
2003). However, since the convergence of the series is generally very slow, the applicability of this 
approach is limited to satellite studies or, in planetary dynamics, to the domain of the so-called 
"hierarchical systems", which are characterized by suitably small values of a. 

In a preceding paper (Michtchenko and Malhotra 2004, hereafter referred to as MM2004), we 
introduced a new semi-analytical approach, consisting of a numerical averaging over the short- 
period perturbations in the mutual interaction of the two planets on coplanar orbits. We obtained 
the planar secular Hamiltonian that allowed us to overcome the limitations of the classical theories. 
Our approach revealed to be suitable for the study of the numerous extra-solar systems whose 
orbital characteristics (e.g. eccentricities and mutual distances) are beyond the range of the classical 
theories. We have investigated the topology of the phase space of the two (coplanar) outer planets, 
c and d, of v Andromedae system (Butler et al. 1997, Butler et al. 1999). Our study permitted 
to explain such common features of planetary systems as the coupled apsidal line motion with 
aligned and anti-aligned pericenters. In addition, the analysis allowed us to discover new dynamical 
features in planetary systems, namely the existence of a nonlinear secular resonance in the very- 
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high-eccentricity region. 

The goal of the present paper is to extend that study to the non-coplanar secular motion of 
the planets. It is known that, at present, the majority of the discovered extra-solar systems are 
spatially unresolved. In fact, the spectroscopic radial velocity technique measures only the line-of- 
sight component of the velocity of the star and the fitting of Keplerian ellipses is unable to detect 
the inclinations and nodes of the planetary orbits. A system must be observed over a large number 
of orbital periods, in order to allow a full N-body fit able to determine all the orbital parameters 
of the system (Ferraz-Mello et al. 2005). The exceptions are the multi-planet systems of strongly 
interacting companions involved in mean-motion resonances. Presently, only the two planets, b 
and c, around PSR B1257+12, which are near the 3:2 mean-motion resonance, have the individual 
orbital inclinations dynamically determined (Konacki and Wolszczan 2003). 

Although the spectroscopic observational data do not constrain the orbital inclinations, dy- 
namical stability considerations can still provide constraints on the individual planetary masses 
and inclinations. For instance, for the v Andromedae system, dynamical arguments suggest that 
the mutual inclination between planetary orbits is less than 20° and that the mean inclination to 
the plane of the sky is greater than 13° (Stepinski et al. 2000, Chiang et al. 2001, Lissauer and 
Rivera 2001). These limits are consistent with the results of the analysis of Hipparcos astrometric 
measurements, yielding an estimate of the true mass of the outermost companion and constraining 
its orbital inclination to the line of sight to be greater than 11° (Mazeh et al. 1999). Also, no 
transits of the innermost planet were found in the photometry observations, limiting its orbital 
inclination to less than 83° (Lissauer 1999). 

It is expected that astrometric and photometric techniques contribute to the task of complete 
characterization of known systems (e.g. the masses and the spatially resolved orbital elements). 
In analogy with the diversity observed in the planar orbital parameters of the known extra-solar 
systems, the elements to be determined may show spatial characteristics that are hard to explain 
within the context of current theoretical models based on the Laplace-Lagrange secular theory 
(Veras and Armitage 2004). We hope that this paper may help to face future challenges. 

The non-linear model developed in this work describes the three-dimensional motion of the two 
planets revolving around the same star. It is based on the same semi-analytical approach used in 
MM2004, consequently, it preserves all advantages of the 2-D model. The model is applied to the 
study of the topology of the 3-D phase space of one planetary system with the dynamical parameters 
of the v Andromedae (planets c and d). The results, presented and discussed in detail in Sections 4 
and 5, are in good accordance with those obtained through direct numerical integration of equations 
of motion. This is mainly due to the fact that the model presented here has no restrictions on the 
magnitude of the planet eccentricities, inclinations and on the planetary spacing. Moreover, the 
present model allows us to characterize the phase space of a given system through the calculation 
of energy levels and families of periodic orbits, independently of high-cost direct integration of 
the exact equations of the motion. The results are completed with direct integrations; it is worth 
emphasizing that the combined use of the two distinct approaches shows very good agreement. 

The paper is organized as follows. In section 2, we introduce the 3-D model and describe the 
basic concepts of the 3-D general three-body problem. They allow us to reduce a number of degrees 
of freedom of the Hamiltonian system to 2. In section 3, we show that the secular behavior of the 
system can be characterized through the geometrical pictures of its phase space. The important 
point is that their construction does not involve the time-consuming numerical integrations. In 
sections 4 and 5, we apply our theory to the outer v Andromedae system, assuming the edge-on 
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observed planets c and d. The topology and dynamical maps of the totally reduced (to 2 degrees 
of freedom) system on the representative subspaces of initial conditions are presented in section 4, 
together with the study of the dependence of the secular dynamics on the magnitude of the Angular 
Momentum Deficit. In section 5, we study the partially reduced (to 3 degrees of freedom) system 
and the dynamical dependence on the initial mutual inclination between the orbital planes of the 
planets. In section 6, we discuss the limits of the applicability of our secular model and explore the 
neighborhood of the v Andromedae system searching for the effects of mean-motion resonances in 
this region. We provide a summary of our results in section 7. Finally, in the Appendix we briefly 
describe the techniques used throughout the paper. 



2 The 3-D model 

The Hamiltonian which describes the motion of the 3-body system in the astrocentric reference 
frame can be written, using Poincare canonical variables, as 

pf mm'j Gm 1 m 2 , (prft) m 

-1 L ^ ' V ✓ 



1= 



~ ~ ^ " direct part indirect part 

Keplerian part 1 

where fi and pi = mi ^ are the planet position vectors relative to the star and their conjugate 
momenta, respectively (pi are the position vectors relative to the center of gravity of the three- 
body system). G is the gravitational constant, /ij = G(mo + rrii), m! i = mjmo/(mo + rrii) and 
A = \ fi — r*2|. Hereafter, the indices % = 1, 2 stand for the inner and outer planets, respectively. 

Associated with the Keplerian part of the Hamiltonian, a set of mass-weighted Delaunay's 
elliptic variables is introduced as 

Mj=mean anomaly, Li=m[yJjLi<ii, 
Wj=argument of pericenter, Gi=Li\Jl — ef, 
r2j=longitude of node, Hi=Li^Jl — ef cos ij, 

where aj, ej and Li are the canonical astrocentric semi-major axes, eccentricities and inclinations 
of the planets, respectively. The relationship between canonical and usual osculating heliocentric 
orbital elements is described in detail in Ferraz-Mello et al. (2004). In the canonical variables, the 
Hamiltonian (1) can be written as 



^/M-rn'j 3 G mi ra 2 



T~t — — % n - 9 R(Li,Gi, Hi, Mi,u)i,Qi), (2) 



where the sum describes Keplerian motions of the planets and R is the disturbing function. 

To study the secular behavior of the system, we perform numerically the averaging of the 
Hamiltonian (2) with respect to the mean anomalies of the planets. The description of the numerical 
averaging procedure can be found in (MM2004). The secular Hamiltonian is then defined by 

H sec = --—2 / mi m2 R{Li, G i} Hi, Mi, ui, ik) dM lC lM 2 , (3) 

(27TJ Z Jo Jo (12 
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where the contribution of the Keplerian part is constant and therefore need not be considered. 

After the elimination of the short periodic terms, the averaged Hamiltonian TC SCC does not 
depend on Mi and M 2 ; consequently, L\ and L 2 (thus a\ and a 2 ) are constant in time. It is 
worth emphasizing that, due to the scissors' kind of the averaging used, the semi-major axes are 
constant only up to the first order in the planetary masses; this consideration is important to bear 
in mind when the results of the model are compared to those obtained through purely numerical 
integrations. 

The secular Hamiltonian H sec given by Eq. 3 has four degrees of freedom and eight variables re- 
main in the averaged problem. However, based on the conservation of the total angular momentum, 
we can reduce the problem by two degrees of freedom. This reduction is known as the elimination 
of nodes (Jacobi 1842). 

For this purpose, we perform the canonical transformation of the nodes to 



3_ ni+n 2 

i\— 2 , 

3 Q1—Q2 



and, accordingly, 



J\=Hi + H 2 =LiJl - e{ cos h + L 2 Jl - e 2 cos I 2 , 



J 2 =H\ - R 2 =L\J\ - e{ cos I\ - L 2 J 1 - e\ cos I 2 . 



Since the averaged disturbing function R depends on the nodes only through AS7 = Q\ — Q 2 
(Brouwer and Clemence 1961), the angle 9\ is cyclic in the expression of R; hence, its conjugate 
action J\ is a constant of motion. 

In terms of canonical astrocentric elliptical elements, the components of the total angular mo- 
mentum C are 



Cx=J2i=i L i\ 1 - e i sin/jsin^j, 



C y =J2i=i L i\J 1 - e i sin Ji cos Hi, 



The direction of the constant total angular momentum defines a plane orthogonal to C, which is 
an invariant of the problem (this plane is known as the Laplace invariable plane). The choice of 
the Laplace plane as the reference plane of coordinates implies that C x = C y = and C z = \ \C\ \ = 
C = const. 

If we assume that both planets orbit around the central star in the same sense, we obtain from 
the condition C x = C y = that the planetary eccentricities and inclinations are such that 

Gisin/i = G 2 sin/ 2 - (4) 

The consequence of this condition is that, for C z > 0, the variation of the inclinations of the planets 
is confined between zero and 90° , with a degeneracy of the problem at Ii = 90° . The longitudes of 
the ascending nodes referred to the invariable plane are coupled in the following way: 

fii = n 2 ± 180°. 

Two immediate consequences of this condition are 1) the system is described by the anti- aligned 
orientation of the nodal lines and 2) the problem is invariant under rotation of the angle Af2 by 
180°. 
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2 n2Mn (5) 



Based on the invariance of the direction of C, we can partially reduce the system by one degree 
of freedom (Malige et al. 2002). In addition, using the invariance of the norm of the total angular 
momentum, we can perform the reduction by one more degree of freedom. Indeed, the condition 
given by C = const constrains the action variables of the problem, in such a way that 

J\=H\ + H2=C, 
J 2 =H l - H 2 =(G( - GD/C. 

Introducing the explicit dependence on C in the Hamiltonian (3), we perform the total reduction 
and obtain, as a result, a two-degrees-of-freedom model. We note that the secular behavior of the 
totally reduced problem parameterized by C is independent upon the mutual inclination between 
the two orbits, I mu t, whose value is defined by the planetary eccentricities, for fixed masses and 
semi-major axes, as 

C'-Gl-Gl 

COS imut = 



2GiG 2 

Instead of the total angular momentum C, the use of the Angular Momentum Deficit (AMD) 
is sometimes more appropriated. This quantity is defined as (Poincare 1897, Laskar 2000) 

2 , 

AMD = L 1 + L 2 -C = ]rL;(l " V 1 " e i cosI i) ( 6 ) 

i=i 

and has a minimum value (zero) for circular co-planar orbits and increases with the eccentricities 
and inclinations. 

In this work, we use both total and partial reductions to study the 3-D dynamics of the planetary 
system. In the first case, we choose a pair of independent action variables and the corresponding 
two angles. The other pair of actions can be easily obtained using Eqs. (5), for a given value of C. 
Here, we opt for two different pairs of the action variables: one is (Gi,G 2 ) (that is, e\ and e 2 ) and 
other is (Gi,H±) (that is, e± and Ii). In the case of partial reduction, the system has three degrees 
of freedom and the set of the variables used consists of the eccentricities of the planetary orbits 
and their mutual inclination. The individual inclinations of the planetary orbits are obtained using 
Eq. (4). 

In the following sections, we will apply the model to a system with the dynamical characteristics 
of the system v Andromedae, restricted to the planets c and d. We adopt the following parameters: 
the masses M = 1.3M , mi sin % = 1.83Mj up and m 2 s'mi = 3.97Mj up , for the central star, inner and 
outer planets, respectively; the semi-major axes, a\ = 0.83 AU and a 2 = 2.56 AU. The description 
of the secular behavior of the system is completed by the initial values of the planet eccentricities, 
ei = 0.252 and e 2 = 0.31 at Azu = 0; formal uncertainties in the eccentricities are of the order of 
20% (Ford 2005). 

Since the observational data do not constrain the individual orbital inclinations, we assume 
that both planetary orbits are observed edge-on, so sini = 1 for each planet, where i is the orbital 
inclination to the plane of the sky. This system will be referred to, along this paper, as the "edge- 
on outer v Andromedae system". The spatial configuration of this system is represented by two 
coplanar planetary orbits placed in the invariable plane, whose inclination i to the sky plane is 90°. 
In this case, the total Angular Momentum of the system is tangent to the sky and, up to second 
order in masses, its modulus can be calculated using the formula: 

2 J 

C = ^mi^JGmo a* (1 - ef). 
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For the fixed masses and semi-major axes of the edge-on outer v Andromedae system, the Angular 
Momentum Deficit (6) is 3 x 10~ 3 , in units of the solar mass, astronomical unit and year. It is 
interesting to note that the edge-on configuration is characterized by a minimal value of the total 
Angular Momentum. Its magnitude increases with the rate 1/sini, due to the increasing planetary 
masses. The AMD also grows in this case. 

We can show that, up to second order in masses, the secular behavior of the systems with 
small I mut is independent on the inclination of the invariable plane to the plane of the sky. This 
is a consequence of the fact that the main features of secular motion of planetary systems are 
independent on the individual values of the planetary masses, but only on their ratio. To show 
this, we first remind that the secular dynamics is defined by such global quantities as total energy 
and the total angular momentum. Then, we point out that, up to second order in the masses, the 
secular energy (Eq. 3) can be divided by the factor Gmim2/a2, while the total Angular Momentum 
Deficit (Eq. 6) can be divided by m2\/Ga2- As a result, we have one problem which depends only 
on the mass ratio and the semi-major axis ratio. 

In the case of coplanar orbits, the mass ratio may be determined from the observations, being 
not affected by the indetermination of the individual masses. As a consequence, the results of the 
edge-on outer v Andromedae system obtained in the study of the coplanar dynamics (MM2004) for 
AMD = 3 x 10 -3 are valid for non edge-on systems with AMD = 3 x 10~ 3 / sini If I mut / 0, the 
determination of the mass ratio is affected by the unknown factor sin i\j sin 22, where i\ and 12 are 
the individual inclinations to the sky plane of the inner and outer orbits, respectively. However, if 
I mut is small, the factor sinii/ sin ^ can still be approximated by 1 that allows us to extend the 
results in the same way as described above. 

Finally, throughout the paper, we will compare the results derived from our semi-analytical 
approach with direct numerical orbit integrations. 

3 Geometrical pictures of the 3-D secular dynamics 

In this section, we present the topology of the phase space of the Hamiltonian system given by 
Eq. (3) and analyze its relationship with planetary dynamic. 

First, to visualize the dynamical features of the fully-reduced system, we introduce a repre- 
sentative plane of initial conditions. The space of initial conditions of the two-degrees-of-freedom 
Hamiltonian system given by Eqs. (2)-(3) is four-dimensional, but the problem can be reduced to 
the systematic study of initial conditions in which two of them are kept constant. This plane may 
be chosen in such a way that all possible configurations of the system are included and thus all 
possible regimes of motion of the system under study can be represented on it. 

We define the angular variables of the representative plane in the following way. The first of them 
comes from the studies of the planar secular problem (MM2004): it is Azu = zu\—zu2 = uji— a>2+Af2, 
where w\ and W2 are the longitudes of the planetary pericenters and AQ = ±180°. We know that 
the angle Azu can either circulate or oscillate about or 180°. In both cases, it goes through either 
or 180° for all initial conditions. Hence, without loss of generality, the angular variable Azu can 
initially be fixed at or 180°. 

The choice of the second angular variable, 2u>\, is based on the particular property of the 
planetary disturbing function that requires the same parity of the indices preceding the planet 
arguments, uj\ and uj2 (Brumberg 1995). Indeed, using this property, we can re-write a generic 
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periodic argument of the Hamiltonian function (3), koj\ + I u>2 + mAQ, as 



— ^— (2u 1 )-l Aw + (m + l) AQ, 

where AO, = ±180° and k, I and m are integers which vary within the interval (— 00, +00). Our 
experimental studies have shown that the angular variable 2 u\ can also either circulate or oscillate 
about 180°, which means that the argument of pericenter of the inner planet, cj±, either circulates 
or oscillates about ±90° (the sign depends on the chosen initial values of angle variables). In this 
case, it always goes through either or ±90° for all initial conditions. Hence, the initial values of 
2u>i can be also fixed at or 180°. 

Now, the dynamics of the twice-reduced Hamiltonian (3) can be represented on the plane of 
the initial values of the pair (ei cos Aw,e2 cos 2u>i), where e± and e2 are the averaged planetary 
eccentricities and the angles Aw and 2o>i are restricted to take the values or 180°. This plane 
will be hereafter referred to as (ei, e^) representative plane. The information on the initial values 
of the angular variables is given by the coordinate signs: positive (negative) values on x-axis 
correspond to Aw = (180°), while positive (negative) values on y-axis correspond to 2uji = 
(180°). 

Figure 1 shows the level curves of the energy given by the averaged secular Hamiltonian (3) on 
the (ei,e2)-plane. The level curves were calculated using AMD = 3x 10~ 3 , which corresponds to the 
edge-on outer v Andromedae system. The thick curve shows the boundary of the energy manifold 
defined by the chosen value of AMD. The distinct spatial orientations of the system describing the 
lower and upper half-planes originate the discontinuity in the representative (ei, ^2) — plane, in the 
transition between these half-planes. 

The geometrical analysis of the secular behavior is based on the fact that the secular energy 
H sec and AMD are both conserved along one solution. Two other facts are also important: (1) 
Independently of whether the motions of Aw and <jj\ are circulations or oscillations, all solution pass 
through the conditions sinAtu = (Aw = and/or 180°) and sin2cji = [oj\ = and/or 90°); 
(2) As shown by numerical simulations, the variables e\ and ei [1\ and I2) reach their maximal and 
minimal values at the condition cos Aw = (cos 2uj\ = 0) , allowing us to estimate the ranges of 
the eccentricity variation for both planets. As a consequence, an individual quasi-periodic solution 
intersects the representative plane at four points and all points must belong to one energy level. Of 
these points, one point is the chosen initial condition, and the other three point are its counterparts. 
The exceptional cases are: a) a stationary solution, which appears as a fixed point on the plane; 
b) periodic solutions, in which one of the angles is fixed, intersect the plane at only two points; 
and c) orbits of chaotic motion, which, due to diffusion processes, intersect an energy level on the 
(ei,e2)-plane at an arbitrary number of points. 

The qualitative behavior of the angular variables, Aw and uj±, can be assessed from the geo- 
metrical analysis of the representative plane in Fig. 1. Concerning the behavior of the angle Aw: 
if all intersections of an orbital path are located at the right-hand side of the (ei,e2)-plane, the 
angle Aw oscillates around zero. Conversely, if all points are located at the left-hand side of the 
plane, Aw oscillates around 180°. Finally, when intersections are found at both half-planes, the 
angle Aw is circulating. 

In what concerns the behavior of the angle uj\, if all points of an orbital path are located at 
the lower half-plane, the angle uj\ oscillates around ±90°. When the intersections are found in 
the lower and upper half-planes, the angle uj\ is circulating. The level of the maximal secular 
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energy, H sec = —1.02370 x 10 , appears in Fig. 1 as a fixed point with coordinates e\ = —0.33 
and e2 = —0.02, corresponding to I mu t = 42°. 9. This is a stationary solution of the secular 
Hamiltonian, which is characterized by constant eccentricities and inclinations. In the close vicinity 
of the stationary solution, the angles Aw and uj\ (consequently, LO2) are oscillating around 180° 
and ±90°, respectively. 

We show in Fig. 1 six levels with decreasing energy denoted by letters from a to f. Several 
examples of intersections of planetary paths with the (ei,e2)-plane are shown by different symbols. 
The level a (H SC c = —1.02377 x 10~ 4 ), close to the stationary solution, is located at the quadrant 
of the plane corresponding to conditions Aw = 2uj\ = 180°. The possible solutions in this case are 
oscillations of both angles Aw and lo\ (and 102) around 180° and ±90°, respectively. Four crosses 
along the level show the intersection of one solution with the representative plane. 

Along the level b, with H scc = —1.0247 x 10~ 4 , the energy level splits into two unconnected 
branches: one is in the lower-left quadrant, while the other is in the lower-right quadrant; both are 
far away from the origin of the plane. The location of this level on the lower half-plane indicates 
that all solutions have u>i oscillating around ±90°. On the other hand, the secular angle Aw 
is allowed now to circulate or oscillate around 180°, depending on the initial conditions. One 
solution on this level, presented by four crosses, has the angle Aw in retrograde circulation. Initial 
conditions leading to solutions with Aw oscillating around 180° with a very small amplitude are 
shown by full circles. As we will see later, the passage from oscillation to circulation of Aw in 
this case is merely kinematical and there is no real separatrix between two modes of motion, which 
evolve continuously from one to another. To make this fact clear in what follows, we will use the 
composite word "circulation/oscillation" to indicate the regime of motion where the two kinds of 
motion co-exist, but are not topologically separated. 

Along the level c, with H scc = —1.0257 x 10~ 4 , most initial conditions are still characterized 
by the libration of uo\ around ±90° and circulation/oscillation of Aw (one example of this kind of 
motion is given by four crosses). The exceptions are solutions in the close vicinity of the origin of the 
plane, where the transition between the oscillating and circulating modes of motion of u\ occurs. 
To better illustrate this feature, we show in Fig. 2 the surface of section and dynamic spectrum 
constructed along the energy level c. The phase space (top panel) is dominated by the regime of 
wi-libration, with the fixed point at eicos2u;i = —0.29, corresponding to e2Cos2u;i = —0.1 and 
^mut = 41°. 8. A separatrix between circulation and libration of u>\ appears when e\ cos2u;i ~ —0.1 
in Fig. 2 and a new regime of motion appears near the origin: this new regime correspond to a 
true secular resonance of the angle Aw, which librates around 180°. The dynamic spectrum in 
Fig. 2 bottom shows its structure with separatrix located at e\ cos2wi = —0.1 (e2 cos2lji = —0.246 
and 7 mut = 34°.7). 

This new dynamical feature is accentuated along the energy level d, with TL sec = —1.027 x 10~ 4 . 
The level is now a continuous curve on the lower half-plane, with a new branch appearing on 
the upper half-plane, where oj\ = 0°. This geometry is also characteristic for the levels e, with 
H scc = -1.029 x 10~ 4 , and f, with H sec = -1.038 x 10~ 4 , in Fig. 1. Along the energy level d, 
we have detected initial conditions leading to the solutions of three distinct regimes of motion, 
depending on the initial conditions. They are 1) a libration of u\ around ±90° and a circulation 
of Aw (shown by crosses); 2) a libration of Aw around 180° and a circulation of uj\ (full circles); 
3) a circulation of both angles (triangles). The regions of transition between domains of different 
regimes are characterized by chaotic motion; one solution of chaotic motion is represented by star 
symbols randomly distributed along the energy level d. To assess the type of motion of one solution, 
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it is necessary to add an information obtained by direct numerical integrations of the equations of 
motion. It will be presented in the next section. 

4 Topology and dynamical maps of the totally reduced system: 
application to the edge-on outer v Andromedae planets 

In this section we investigate the whole phase space of the Hamiltonian system given by Eq. 3, with 
fixed values of AMD. Our goal is to identify all possible regimes of its secular motion. For this 
task, we apply the model to one system whose masses and semi-major axes are those of the edge-on 
outer v Andromedae system. The results are presented in the form of topology and dynamical 
maps. The topology maps were calculated using the semi-analytical approach described in the 
previous sections, while the dynamical maps were obtained through numerical integrations of the 
exact equations of motion. 

To fully represent the secular dynamics, we have chosen two subspaces of initial conditions: one 
is the (ei,e2)-plane introduced in the previous section, and the other is a similar plane on which 
e<i is substituted by the inclination of the inner planet over the invariable plane, the (ei,ii)-plane. 
In both cases, the initial values of the angular variables of the problem, Aw and 2u>i, were fixed at 
either zero or 180°. The information on the initial values of the angular variables is given by the 
coordinate sign, as described in the previous section. 

It is easy to show that both phase planes are equivalent: the solutions on one of these planes 
may be obtained from those on the other through Eqs. (5), for a given AMD. Nevertheless, we 
present two planes due to the fact that the (ei,e2)-plane is appropriate to display the details in high 
planetary inclinations, but the low inclination behavior is confined to the very narrow vicinity of the 
borders. At variance, the (ei,ii)-plane shows clearly the details of the low inclination dynamics. 

4.1 Phase portrait of the edge-on outer v Andromedae system 

Figure 3 shows the topology (left) and dynamical maps (right) on the (ei,e2)~plane (top) and (e\,I\)- 
plane (bottom), both calculated for AMD = 3 x 1CT 3 , which corresponds to the edge-on outer v 
Andromedae system. The location of the coplanar planets c and d is shown by a star symbol on all 
panels. The level curves of constant energy are plotted with solid lines in the topological maps. The 
mutual inclination between the initial orbital planes is given by dashed lines on the (ei,e2)-plane. 
On the (ei,ii)-plane, the dashed lines show the levels of e2 = const. For AMD = 3 x 10~ 3 , the 
system reaches the maximal value of the mutual inclination at 47°. 4, when e\ = e2 = 0, while the 
maximal value of the eccentricity of the outer planet is 0.365, at e\ = I\ = 0. 

The main families of periodic orbits are plotted in the topological maps of Fig. 3 (left) by color 
dots: The red dots are associated with solutions of the Hamiltonian (3) in which the angle Aw 
is oscillating, while the blue dots to solutions with librating u\. The stability of the periodic 
solutions was inferred from the features of the dynamical map (see description of the dynamical 
maps below). The character of motion on each orbit is assessed by the identification of the regimes 
of motion of the system by means of numerical integrations. Note that the oscillation of the 
secular angle Aw does not imply necessarily a resonant behavior, but may be just a kinematical 
continuation of the circulation regime of motion (see MM2004). It is worth mentioning that, except 
these circulation/oscillation solutions, the families present loci (stable by large dots and unstable 
by small dots) of the secular resonances in the phase space of the system and are in very good 
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agreement with the location of the secular resonances obtained by numerical integrations. The 
calculation of the periodic solutions provides an important information about a dynamical system, 
without time-expensive numerical integrations. 

In Fig. 4 we show four examples of periodic orbits: two orbits correspond to the red curves 
in Fig. 3 and are, respectively, stable (left-top) and unstable (left- bottom). In the stable solution, 
the angle Aw librates around 180°. The other two orbits correspond to the blue curves and are, 
respectively, stable (right-top) and unstable (right- bottom). In the stable solution, the angle uj\ 
librates around ±90°. 

In the construction of the dynamical maps of Fig. 3 (right), a grid of 72 x 79 initial conditions 
was defined on the (ei,e2)-plane, with steps Aei = 0.02 and Ae2 = 0.01. The initial inclinations 
were obtained from Eqs. (5), for AMD = 3 x 10~ 3 . The initial value of the angular variables used 
were: from the set I in the right-top quadrant, the set II in the left-top quadrant, the set III in 
the left-bottom quadrant and the set IV in the right-bottom quadrant of the (ei,e2)-plane (see the 
definition of the sets in the Appendix). On the (ei,/i)— plane, a grid of 72 x 81 initial eccentricity 
and inclination of the inner planet was defined with steps Aei = 0.02 and AI± = 1°. The initial 
eccentricity and inclination of the outer planet were obtained using Eqs. (5). The initial angular 
orbital elements of the planets are the same as used in the previous case. 

The shading scale used in the dynamical maps in Fig. 3 (right) is related to the degree of 
stochasticity of the solutions: the lighter regions in the dynamical maps correspond to initial 
conditions of regular motion, darker tones indicate increasingly chaotic motion. The limit of the 
domains of allowed and forbidden motion can be calculated using Eqs. (5): with conditions I\ = 
I2 = on the (ei,e2)-plane and e2 = h = on the (ei,Ii)-plane. Outside, there are no solutions 
of the Hamiltonian system (3), for AMD = 3 x 10~ 3 . 

The analysis of the structure of the dynamical maps reveals interesting topological properties 
and reflect important features of the 3-D secular dynamics, when the orbital elements of the edge- 
on outer v Andromedae system are used. The (ei,e2)~ and (e±,Ii) -planes are dominated by a gray 
background of regular secular motion of the system. The narrow white strips coincide with the 
location of the stable periodic solutions of the secular system. The domains of chaotic motion 
appearing as black regions in Fig. 3 (right) are associated with separatrices of different regimes of 
secular motion (and unstable periodic orbits). 

The careful analysis of the results of the numerical simulations, always accompanied by the 
study of the geometry of the energy levels and periodic solutions, allowed us to identify various 
different regimes of motion of the system. Their domains are marked by 1-4 on the dynamical 
maps in Fig. 3 (right). The domains 1 are regions of motion characterized by the coupled variation 
of the eccentricity and inclination of the inner planet and the libration of the angle uj\ around ±90°. 
They are located on the lower half-planes at high mutual inclination (above 35°). The white strips 
follow the location of the periodic resonant solution in these regions. This regime of motion is often 
referred to as Kozai resonance in the literature. In this work, we designate this resonance as e—I 
coupling, or Lidov-Kozai resonance, to pay homage to the Russian scientist, who first discovered 
this dynamical phenomenon (Lidov 1961). The stationary solution of the secular Hamiltonian (3), 
together with the paths of energy levels from a to c shown in Fig. 1, belong to the domains of this 
resonance. Inside the Lidov-Kozai resonance, the angle Aw shows either a retrograde circulation 
or an oscillation around 180°. Two examples of the periodic orbits (stable and close-to-unstable) 
in this regime of motion were shown in Fig. 4 (right ). Note that the unstable orbit escapes from 
the Lidov-Kozai resonance after a while. 
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The domains of the Lidov-Kozai resonance are delimited in the phase space by regions of strong 
chaotic motion, which appear as black zones in the dynamical maps. In their close vicinity, at lower 
mutual inclinations, a new resonant regime of motion appears; its domains are labeled by 2 in Fig. 3. 
This resonance is characterized by the true libration of the secular angle Aw around 180° and the 
prograde circulation of both angles oj\ and oj 2 - The family of the stable periodic solution of this 
resonance can be observed as the narrow white strips inside the domains 2. Two examples of the 
periodic orbits (stable and close-to-unstable) in this regime of motion were shown in Fig. 4 (left 
column). 

The regime of motion in the region 2 is better illustrated on a surface of section. In Fig. 5, 
we show the surface of section and dynamic spectrum constructed along the energy level TC sec = 
-1.033 x 10" 4 . The dominating regime of motion at this energy is the secular resonance of Aw 
described above, with the stable fixed point at e\ cos Aw = —0.385 [ei cos 2oj\ = —0.187 and I mut = 
32°. 7). The unstable fixed point of this resonance is located at eicosAtu = 0.433 (e2cos2oji = 
—0.147 and the same mutual inclination). 

We can see in Fig. 5 that the separatrix involving the secular resonance divides the whole domain 
of motion into two zones: the inner zone, close to the origin, and an outer zone. The inner zone 
shows a complex dynamical structure characteristic of the presence of secondary resonances. The 
dynamic spectrum of the solutions with initial conditions e\ sin Aw = (Fig. 5 bottom) reveals the 
existence of islands of regular motion inside the sea of chaos in this region. The secular angles inside 
the inner zone are in circulation: retrograde for Aw and prograde for u\ and The domains of 
this regime of motion are marked with the label 3 in the dynamical maps of Fig. 3. 

In the outer zone, the circulation of Aw inverts its direction. This regime of motion is charac- 
terized by prograde circulation of all angles: Aw, oj\ and uji- The domains of this regime of motion 
are marked with the label 4 in Fig. 3. This regime of motion extends to mutual inclinations close 
to and is known from the study of dynamics of the planar planetary model (MM2004). For some 
initial conditions, the angle Aw oscillates around or 180°; the regions of oscillation of Aw in the 
domains 4 are seen as the white strips on the dynamical maps in Fig. 3. It should be emphasized 
that the difference between circulation and oscillation of Aw in this region is merely kinematical 
and there are no true separatrices between the two modes of motion (MM2004). In other words, the 
apsidal alignment (or anti-alignment) in this case is a simple circulation around a center displaced 
from the origin and not a resonant motion. The 3-D dynamics of the system in the domains 4 is 
regular and nearly similar to that of the coplanar case, even for mutual inclinations as large as 30°. 
However, darker tones on dynamical maps in Fig. 3, also appear indicating increasingly chaotic 
motion. The explanation for this feature is the proximity of the system to mean-motion resonances 
of high order (see Section 7). 

Finally, the regimes of motion appearing in Fig. 3 are summarized in Table I (the regime 5 will 
be discussed in the next section). 

4.2 Dependence on the Total Angular Momentum Deficit 

The secular behavior of the fully-reduced system is determined by the adopted value of the Total 
Angular Momentum Deficit. In the previous section we have shown the main features of the 
dynamics obtained for AMD = 3 x 10~ 3 , which corresponds to the edge-on outer v Andromedae 
system. However, the exact value of this quantity is unknown, mainly, due to the fact that the 
observational data contain no information on the mutual inclination between the planets and the 
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inclination of the invariable plane. For this reason, in this section, we investigate the evolution of 
the secular dynamics function of AMD. 

We have chosen two distinct values: one is smaller (AMD = 1 x 10 -3 ), and other is larger 
(AMD = 8 x 10~ 3 ), than AMD of the edge-on v Andromedae system, notwithstanding that AMD 
of the actual outer v Andromedae system cannot be smaller than the edge-on value. The smaller 
value of AMD is however considered in this work to complete the dynamical picture of the secular 
system. At variance, larger values of AMD are not impossible. For instance, it is enough to assume 
that both orbits are edge-on (to keep the same masses), but have a non-zero mutual inclination. 
(The chosen value is purposely large and corresponds to I mut « 60°.) The main results obtained 
are given in the following. 

For small values of AMD, the dynamics of the system is very similar to that of the planar 
case. This is due to the fact that the values of the mutual inclination between the planetary orbits, 
allowed from Eqs.(5), in this case, are small. In Fig. 6, we present the topological picture of the 
phase space of the system with AMD = 1 x 10~ 3 , plotting the energy levels on the (ei,e2)- and 
(ei,ii)-planes. The maximal value of the mutual inclination between orbital planes is 27°. The 
typical signatures of the high-inclination behavior, such as the existence of stationary solutions and 
bifurcations of the energy levels, shown in the previous sections, are not observed in Fig. 6. Only 
periodic solutions corresponding to the secular angle Aw exist. These periodic solutions show the 
centers around which Aw oscillates (see MM2004). 

At variance, for the large values of the Angular Momentum Deficit, the domain of possible 
motion of the system is extended to very high eccentricities and inclinations. Figure 7 shows the 
topology (left) and dynamical maps (right) on the (ei,e2)— plane (top) and (ei,/i)-plane (bottom), 
constructed for AMD = 8 x 10~ 3 . The notations used in this figure are analogous to those used in 
Fig. 3. The system reaches the maximal value of the mutual inclination at 79°. 7, when e\ = e2 = 0; 
the maximal value of the eccentricity of the outer planet is 0.58, when e\ = I\ = 0. The location 
of the fixed point (stationary solution) is at e\ = —0.81 and C2 = —0.1, which correspond to 
imut = 60°. 9. 

The geometry of the representative planes is similar to that obtained for AMD = 3 x 10~ 3 . The 
structure appearing in Fig. 7 left bottom, in the low inclination region along the negative x-axis is 
an artifact. It only exists in the averaged model, disappearing when the full exact equations are 
considered. 

Previous studies of the planar problem (MM2004) have shown that the large values of the 
Angular Momentum Deficit are characterized by the presence of a secular resonance with the angle 
Aw librating around 0. Indeed, the domain of this resonance, denoted with the label 5, is clearly 
seen on the dynamical map of the (ei,7i)-plane in Fig. 7 right bottom; on the (ei,e2)-plane, this 
regime of motion occupies a narrow domain near the border. 

The high- inclination regimes of motion 1 and 2, discussed in the case of Fig. 3, are also present 
in Fig. 7. The regime 1, that is, the Lidov-Kozai resonance, dominates the lower half-planes. The 
location of the periodic orbits of this resonance is in good agreement with that calculated from 
our 3-D model (blue curves on the topological maps). The regime 2, which is the high- inclination 
secular resonances of Aw, dominates the upper half-planes: one in the right side, where Aw librates 
around 0, and the other in the left side, where Aw librates around 180°. 
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5 Partially reduced dynamical system: dependence on the mutual 
inclinations 



In the previous sections, the analysis of the 3-D dynamics was done with a system reduced to 
2 degrees of freedom, whose evolution is constrained by a fixed value of the Angular Momentum 
Deficit. In this section, we consider the system only partially reduced to 3 degrees of freedom, where 
the direction of the angular momentum is taken into account, but not its modulus. In this case, 
we vary three orbital parameters of the problem; they are chosen as both planetary eccentricities 
and the mutual inclination of the planetary orbits. The study is done by means of dynamical maps 
and characteristic curves of periodic orbits of the Hamiltonian system given by Eq. 3. 

Figures 8 and 9 show dynamical maps on the (ei,e2)-plane, obtained for three different values 
of the initial mutual inclination, namely 5°, 15° and 45°. In the construction of the maps, a grid of 
100 x 61 initial eccentricities was defined on the (ei,e2)-plane. The angular variables used were: the 
set III, in the left-hand side of the panels, and the set IV, in the right-hand side of the (ei,e2)-plane 
(see the Appendix). 

For I mut = 5° and 15°, the domains of regular motion are characterized by the circula- 
tion/oscillation regime of motion of Aw. The white regions of the phase space correspond to 
oscillation of Aw: either around or around 180°. The regions coded in gray are regions of regular 
motion with circulating Aw. Finally, the regions of highly nonharmonic and chaotic motion (dark 
zones) are domains of initial conditions leading to close approaches of the two planets. 

The hatched regions are regions of large-scale instability followed by disruption of the system 
within the time-interval of 530,000 years. The location of the edge-on outer v Andromedae system 
is shown by a star on both panels in Fig. 8. For these parameters, the system is within the domain 
where Aw oscillates around zero, as has been noted in previous studies (Malhotra 2002, Chiang & 
Murray 2002, Ford et al. 2005). The effect of the mean- motion resonances of high orders may be 
noted as a weak nonharmonic gray feature in the close proximity of the system on both panels in 
Fig. 8. 

The domain of the true secular resonance reported by MM2004 is seen at the high eccentricity 
region of regular motion near the right border. The resonant orbits are protected from close 
approaches by coupled variation of the planet eccentricities and libration of Aw. It is interesting 
to observe that the orbit of the inner planet inside this resonance can reach eccentricities as high 
as 0.95, but the system of two planets continues to be stable (at least, over the time-interval of 
530,000 years). 

The differences between the dynamical features revealed by two planes in Fig. 8 (obtained for 
-fmut = 5° and 15°) are negligible. We conclude that the dynamics of the 3-D system varies slowly 
in the low-inclination region and is similar to the dynamics of the planar system studied in detail 
in (MM2004). This conclusion is in agreement with the results obtained by Lee and Peale (2003) 
for hierarchical planetary systems. 

Significant changes in the dynamics of the system take place for initial mutual inclinations 
larger than 30°. The dynamical map obtained with I nmt = 45° is shown in Fig. 9. The Lidov- 
Kozai resonance is now dominating over the whole domain of regular motion. The analytically 
calculated characteristic curves of the periodic orbits of this resonance (in blue color) are in very 
good agreement with the features of the dynamical map (white strips) . 

The secular angle Aw circulates for almost all initial conditions, but, depending on the initial 
conditions, there are also oscillations around 180°. These oscillations, which are just a kinematical 
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continuation of circulations, occur in the low e 2 region, where the red curve shows their location 
in Fig. 9. In the very high eccentricity region, the red curve indicates a narrow domain of stable 
motion: this is the domain of the true secular resonance characterized by libration of the angle Azu 
around 180°. 

Finally, domains of chaotic motion appear in Fig. 9 as inclined black strips of variable width; 
they are associated with the effects of the mean-motion resonances in the neighborhood of the 
system. 

6 Mean-motion resonances in the neighborhood of the v Andromedae 
Planetary System 

The applicability of the model of secular dynamics presented in this work is restricted to domains 
free from the influence of mean-motion resonances. The averaging procedure of Eq. 3 removes, 
together with short-period terms, all effects of mean-motion commensurabilities; as a consequence, 
the semi-analytical approach can not provide any information about mean-motion resonances and 
dynamical instabilities associated to them. The validity of our secular model, may be analyzed 
from the dynamical maps obtained by numerical integration of the exact equations of the motion, 
where the effects of mean-motion resonances in the neighborhood of the system under study are 
always present. For this reason, we decided to investigate the neighborhood of the v Andromedae 
system looking for the peculiar features of mean-motion resonances. 

Special caution was taken with the sini indetermination in the masses of the planets. At 
variance with the secular case, there are no evidences that the resonant behavior depends only on 
the mass ratio; therefore the dynamical structure of the phase space may be strongly sensible to 
the inclination of the orbits to the sky plane. For this reason, we use two sets of the planetary 
masses: The current one corresponding to the edge-on outer v Andromedae system, with the masses 
mi = l-83Mj up and m 2 = 3.97Mj up , and a second set adopting sini = 0.5, which corresponds to 
the inclination of 30° on the sky plane, with the masses m\ = 2 x 1.83Mj up and ni2 = 2 x 3.97Mj up . 
(Note that the mass ratio is the same in both sets). 

In both cases, the initial orbital elements in the simulations were uniformly distributed on 
the (a2,e2)-plane of the semi-major axis and eccentricity of the outer planet, within the domains 
2.4 AU< a 2 < 2.7 AU (Aa 2 = 0.005 AU) and < e 2 < 0.6 (Ae 2 = 0.01), respectively. The initial 
inclination of the outer planet to the invariable plane was fixed at 1°. The semi-major axis and 
eccentricity of the inner planet were fixed at their current values, a\ = 0.83 AU and e\ = 0.252. 
We have chosen the set I of initial values of angular variables, and, to have the motion referenced 
to the invariable plane, we have calculated the inclination of the inner planet using Eq. (4). 

Figure 10 shows the dynamical maps of the neighborhood of the current position of the outer 
planet d: on the left panel, for the first mass set, and on the right panel for the second set. Once 
again, the shading scale was used to distinguish between the regions of regular and chaotic motion. 
In both panels in Fig. 10, the hatched regions are characterized by large-scale instabilities followed 
by disruption of the system within the time-interval of 530,000 years. The domains of chaotic 
motion are associated with mean- motion resonances. We note the dominating presence of the 5:1 
resonance: this is the resonance of the lowest order (4) in this region. Weaker resonances of higher 
orders appear as black strips of variable width: they are the 11:2, 16:3 and 21:4 resonances of orders 
9, 13 and 17, respectively. The location of the v Andromedae system is marked by a star symbol 
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on both panels. 

We note that the edge-on system shown in Fig. 10 left is far away from the region of the strong 
5:1 resonance; this fact may validate the use of the secular model in the study of this system. On 
the other hand, it is close to the weak 16:3 mean-motion resonance. We believe that the situation 
is not dramatic in this case, although the resonant effects can be still seen on dynamical maps. Due 
to its high order, the 16:3 resonance is very narrow; moreover, due to constraints imposed by the 
conservation of the total angular momentum, this system may survive for a time comparable with 
the age of the star (Michtchenko and Ferraz-Mello 2001). 

The situation may be considered as dramatic in the case of the second set of the adopted 
planetary masses, corresponding to a 30°-inclined invariable plane (Fig. 10 right). In this case, 
the region of large-scale instabilities increases and all mean-motion resonances are shifted towards 
the larger semi-major axes. The system under study is located, now, very close to the strong 5:1 
mean-motion resonance. No doubt that, in this case, the application of the secular model to v 
Andromedae would be hazardous. 



7 Conclusions 

This work presents an extension of the planar secular semi-analytical model introduced in (MM2004) 
to the study of the three-dimensional dynamics of planetary systems. The basic technique is the 
numerical averaging over short-periods of the mutual interaction of the two planets. We emphasize 
that, in the present work, we do not use expansions of the disturbing function in power series of 
the orbital elements nor in Fourier series of the angular variables, such as done classically (e.g. 
Brouwer and Clemence 1961, Murray and Dermott 1999). This means that the secular motion 
of the planets is described very precisely, without any restriction about the magnitude of their 
eccentricities, inclinations and the mutual distance. The only condition for the applicability of the 
model is that the system must be located far from a strong mean-motion resonance. 

Due to the invariance of the total angular momentum in the general three-body problem, the 
3-D secular system is reducible to a two-of-degrees-of-freedom dynamical system (total or Jacobi 
reduction) . The behavior of the totally reduced system is governed by the adopted value of the total 
Angular Momentum Deficit (AMD). Recall that the AMD is equal to zero for circular co-planar 
orbits and increases with increasing eccentricities and/or inclinations. We have shown that, up 
to second order in masses, the 3-D phase space structure of the secular system does not depend 
on the individual values of the planetary masses and semi-major axes, but only upon the ratios 
of these quantities. For a given AMD, the mutual inclination between the planetary orbits is a 
function of these parameters and the planetary eccentricities. For small mutual inclinations, the 
secular dynamics is independent of the orbital inclinations to the plane of the sky. This allows us 
to overcome partially the problem of the mass indetermination. 

Using the non-linear secular 3-D model, we have constructed geometrical pictures of the secular 
phase space of the two-planet system in terms of eccentricities and inclinations. Owing to the 
symmetries in the secular Hamiltonian, the phase space structure can be visualized in representative 
planes of initial conditions, with the initial angular elements Aw and 2uj\ fixed at either or 180°. 
The analysis of the topology of the phase space of the Hamiltonian system given by Eq. (3) allow us 
to estimate the range of the eccentricity /inclination variations without time-expensive numerical 
integrations of the equations of motion. 

We have investigated the whole phase space of the Hamiltonian system given by Eq. 3 looking 
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for all possible regimes of secular motion. For this task, we applied our approach to the specific 
case of the system with a dynamical parameters of two outer planets, c and d, of the edge-on v 
Andromedae. The qualitative study was always supplemented by direct numerical integrations over 
a wide range of initial conditions and the interpretation of the results was always founded on the 
analytical approach. 

The topology of the phase space of the system was investigated by means of several techniques, 
nominally: energy level maps , characteristic curves of families of periodic orbits , surfaces of 
section, dynamic spectra and dynamical maps. 

In the following we summarize the important features of the 3-D dynamics of the secular system. 

1. The low-to- moderate eccentricity and mutual inclination regime of motion (domain 4 with 
e\ < 0.6 and I mu t < 30°). This is a general non-resonant regime, similar to that found in the planar 
case. The systems always exhibit two main modes of secular motion, characterized by circulation 
of Aw or its oscillation around or 180°. There are no real separatrices between the two modes 
of motion, whose solutions evolve continuously from one type to another; the differences between 
these solutions are merely kinematical. The apsidal alignment or anti-alignment in this case is a 
simple circulation around a center displaced from the origin. The arguments of pericenter are in 
regular direct circulation. The edge-on outer v Andromedae system is in this regime of motion and 
its secular behavior is stable; weak instabilities seem to occur only in the vicinity of the high-order 
16:3 mean-motion resonance. 

2. The high eccentricity and low-to-moderate inclination regime (domain 5 with e\ > 0.6 and 
-fmut < 10°) is characterized by large-scale instabilities, due to close approaches of the planets, 
followed by disruption of the system within a few thousands of years. The only surviving solutions 
in this region are those inside the nonlinear secular resonance and are bounded by a zero-frequency 
separatrix. The secular angle Aw librates around and the variation of e\ and e2 is strongly 
coupled. This feature of the secular dynamics of two-planet systems is the same that appeared in 
the studies of the planar problem (MM2004). 

3. The high inclination regimes of motion (domains 1-3 with J mut > 30°): complex dynam- 
ical behavior with the presence of several regimes of resonant motion. The dominating behavior 
is the e-I coupling, or the Lidov-Kozai resonance, characterized by the coupled variation of the 
eccentricity and inclination of the inner planet and the libration of the angle u\ around ±90°. At 
variance with the analogous phenomenon in restricted problems, the variation of the planet incli- 
nations and eccentricities is constrained by the total angular momentum conservation. The large 
eccentricity /inclination excursions induced by the Lidov-Kozai resonance in restricted problems 
can not occur in the planetary problem. A regime of motion with Aw in the secular resonance also 
exists in the high-inclination region. In this case, the secular angle Aw librates either around or 
180°. This regime is new and has no relation with the true secular resonance found in the planar 
problem by Michtchenko and Malhotra (MM2004) 

Finally, the limits of applicability of the model were assessed by the construction of dynamical 
maps. We remind that the validity of the results derived from the secular approach may be 
compromised by effects of mean- motion resonances. We detect the presence of the several mean- 
motion resonances in the neighborhood of the v Andromedae system: one is the low-order 5:1 
resonance, and the others are all of very-high order with no significant effects on the secular motion 
of the system. The system with the masses of the edge-on outer v Andromedae system is located 
far away from the strong 5:1 resonance. 
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8 Appendix. Applied techniques 



In this section we describe the main numerical techniques used in the study of the 3-D dynamics 
of the planetary system. 

Energy level maps. The topology of the phase space of the Hamiltonian system can be studied 
by plotting the energy level curves on representative planes of initial conditions. For this purpose, 
the equation H scc — H* = is solved numerically, for a given value of the energy H* , using a numeri- 
cal procedure for root finding. The secular Hamiltonian function 7i scc is given by Eq. (3). Applying 
the total reduction, we choose a pair of action variables and the other pair is obtained through 
Eqs. (5). The pair of independent variables may be, for instance, the pair (Gi,Gi) (equivalently, 
ei,e2). In this case, for given 7i* and AMD, e\ can be easily obtained as a function of e<i- Varying 
H*, we calculate a family of the energy levels and plot it on the representative (ei,e2)-plane. 

Characteristic curves of families of the periodic orbits. The location of the zero precession 
rates of the angular variables were obtained through the conditions: 

Ul ~ dGi ~ u ' W2 ~~ dG 2 ~ U ' ^ > 

By definition, the difference of the planet periastron longitudes is Aw = u>\ — U2 + AfJ, where 
A£l = 180°. Consequently, we obtain the condition Am = lo\ — C02 = for the zero precession rate 
of the secular angle Aw. The derivatives can be computed numerically for any given point of the 
phase space using the second-order differentiation scheme (Press at al. 1986). 

It should be emphasized that, for the system with two or more degrees of freedom, the conditions 
(7) do not define families of periodic orbits. In this work, we present the curves corresponding to 
tb± = and Aw = 0. They show the domains of possible libration or oscillation of the corresponding 
angles. 

A two-degrees-of-freedom dynamical system is characterized by two fundamental frequencies, 
consequently, a general (quasi-periodic) solution of this system is a composition of two independent 
modes of motion. In this case, the condition of periodicity will be satisfied when the amplitude of 
one of two modes tends to zero. The precise location of the zero-amplitude periodic trajectories 
can be calculated numerically from the position of the fixed points on the surfaces of section. 

Numerical integrations. In order to obtain planetary motions, the exact equations of planetary 
motion were numerically integrated in the framework of the three-body general problem. The accu- 
rate RA15 integrator (Everhart 1985) was used. To compare the results of numerical investigations 
with the results given by the developed model, we use the canonical astrocentric orbital elements of 
the planets, instead the usual osculating astrocentric elements. The transformation between these 
two sets can be found in detail in (Ferraz-Mello et al. 2004). 

The initial angular orbital elements of the inner (i = 1) and outer (i = 2) planets used in the 
dynamical maps were: 

Set I (cosAtu = 1 and cos2uji = —1): 

^=-90° Mi=0° 
lo 2 = 90° M 2 =0° 



20 



Set II (cos Act = — 1 and cos2u;i = — 1): 

cji=90° Mi=180° 
cj 2 =90° M 2 = 0° 

Set III (cos Act = 1 and cos2a;i = 1): 

a;i=-180 o M 1= 0° 
u 2 = 0° M 2 =0° 

Set IV (cos Aw = —1 and cos2wi = 1): 

Ul =0° Mi=180° 
uj 2 =0 o M 2 = 0° 

For all sets, the initial value of AO = Q\ — Q 2 was fixed at 180°. 

In order to apply the theory developed in this paper to the outer v Andromedae system, 
we have adopted the following parameters: the masses m sta r = l-3-Ms un , mi = 1.83Mj up and 
m 2 = 3.97Mj up ; the semi-major axes Gq = 0.83 AU and a 2 = 2.56 AU; and the planet eccentricities 
ei = 0.25 and e 2 = 0.34 at Aw = and 2wi = 0. 

To remove the short-period oscillations (those of the order of the planetary orbital periods), a 
low-pass filtering procedure was implemented on-line with the numerical integration as described in 
detail by Michtchenko and Ferraz-Mello (1995). The numerical integrations were performed over a 
time interval of 524,544 years, which was large enough to allow an accurate and efficient averaging 
of the long-period effects, and also to detect the occurrence of secular resonances. 

Dynamical maps. The orbits of the planets obtained trough direct numerical integrations were 
Fourier-transformed using the standard FFT algorithm. The information contained in the power 
spectra of the orbital elements was used in the construction of dynamical maps (see Michtchenko 
et al. 2002). 

The mapped quantity is the spectral number TV defined as the number of peaks in the power 
spectrum of one calculated planetary orbit above an arbitrarily defined "noise". In this work, we 
consider as "noise" those peaks with amplitudes smaller than 5% of the amplitude of the largest 
peak. The spectral number N is used to qualify the chaoticity of planetary motion in the following 
way: small values of N correspond to regular motion, while the large values indicate the onset of 
chaos. 

Once the spectral numbers N are determined for all the initial conditions on the grid, we plot 
them on the representative planes using a shading scale. The calculated values of N, in the range 
from 1 to 80, are coded by a gray level scale that varied linearly from white (N = 1) to black 
(N = 80). Since large values of N indicate the onset of chaos, the shading scale is related to the 
degree of stochasticity of the initial conditions: lighter regions on the dynamical maps correspond 
to regular motion, darker tones indicate increasingly chaotic motion. 

Surfaces of section and dynamic power spectra. The structure of the phase space of the 
two-degrees-of-freedom system was studied using the standard technique of the construction of 
surfaces of section. 
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Two sections have been chosen for presentation of the totally reduced system: The first is a 
section by the plane sin 2u\ = and its coordinates are e\ cos Aw and e\ sin Aw. The second is a 
section by the plane sin Aw = and its coordinates are e± cos2w! and e\ sin 20^. 

It should be emphasized that, in the construction of the surfaces of section, we have used as 
input the planetary solutions obtained by numerical integration (the application of on-line filtering 
procedure is crucial in this task). For this reason, the erratic scatter of points due to the loss of 
numerical accuracy can be sometimes observed. 

The dynamic power spectrum is a technique which is complementary to the surfaces of section 
and is very efficient to identify such phenomena as bifurcation and chaoticity. A dynamic spectrum 
presents, for each value of a given dynamical parameter (the abscissa of the plot), the evolution of 
the main oscillation modes of the planetary motion as a function of the parameter. Over the domains 
of regular motion, the proper frequencies vary continuously when the parameter is gradually varied. 
When the region of chaotic motion is approached, the frequency evolution shows a discontinuity, 
characterized by the erratic scatter of values. 
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Table 1: The main regimes of secular motion. The behavior of the angles u\ and Aw is identified 
by symbols: L - libration, CP and CR - prograde and retrograde circulation, respectively, O - 
oscillation. 



Label 




Aw 


^mut 


ei 


1 


L (±90°) 


CR / (180°) 


> 35° 


all 


2 


CP 


L (0 and 180°) 


> 30° 


all 


3 


CP 


CR 


~ 30° 


all 


4 


CP 


CP / (0 and 180°) 


< 30° 


< 0.6 


5 


CP 


L(0) 


< 10° 


> 0.6 
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Figure Captions 



Figure 1. Schematic details of the phase space on the representative (ei,e2)-plane of initial condi- 
tions obtained for AMD = 3 x 10~ 3 . Positive (negative) values of e\ correspond to Aw equal to 
(180°); positive (negative) values of e2 correspond to 2cji equal to (180°). The distinct spatial 
orientations of the system describing the lower and upper half-planes originate the discontinuity 
in the transition between these half-planes. Six levels, from a to f, with decreasing energy, are 
shown. The fixed point is a stationary solution of the Hamiltonian (3). The intersections of some 
numerically calculated orbits with the plane are shown by crosses, full circles, triangles and stars 
symbols (see text for further details). 

Figure 2. (Top:) Surface of section defined by sinAro = and AMD = 3 x 1CT 3 . (Bottom): 
Dynamic spectrum corresponding to the solutions whose initial conditions are on the axis e\ = 
of the surface of section. 

Figure 3. Left: In black lines, energy levels of the secular Hamiltonian given by Eq. (3) on the 
(ei,e2)— plane (top) and (ei,ii)-plane (bottom), for AMD = 3x 10 -3 . The signs + or — , preceding the 
variable e±, indicate that the initial values of Aw are zero or 180°, respectively. The signs + or — , 
preceding the variables &2 and l\ indicate that the initial values of 2cji are zero or 180°, respectively. 
In dashed lines, level curves of the mutual inclination on the (ei,e2)-plane and the eccentricity of 
the outer planet (ei,7i) -plane. Periodic solutions of the secular Hamiltonian, corresponding to 
oscillating Aw or librating uj\ are shown by red and blue dots (large for stable and small for 
unstable), respectively. The location of the edge-on outer v Andromedae system is indicated by 
a star. Right: Dynamical maps. The domains of the different regimes of motion are: 1 - the 
Lidov-Kozai resonance, where the angle lv\ is in libration around ±90° and Aw is in retrograde 
circulation; 2 - the secular resonance, where the angle Aw is in libration around 180° and u\ is 
in circulation; 3 - the angles u\ and Aw are in direct and retrograde circulation, respectively; 4 
- Aw is in the general regime of circulation/oscillation of the planar case and both arguments of 
pericenter are in direct circulation. The domains of forbidden motion, with no solutions of the 
Hamiltonian system for AMD = 3 x 1CT 3 , are hatched. 

Figure 4. Four of the periodic solutions shown in Fig. 3. Each panel presents two curves: thick 
curve is defined on the (e\ cos Aw x e\ sin Aro)-plane and thin on the (I\ cos 2uj\ x I\ sin 2wi)-plane. 
Left column: The stable (top) and unstable (bottom) orbits with the secular angle Aw librating 
around 180° and/or circulating. Right column: The stable (top) and unstable (bottom) orbits with 
the angle uj\ librating around ±90° and/or circulating. 

Figure 5. Same as Fig. 2, but with the condition sincji = 0. 

Figure 6. Same as Fig. 3 (left column), but for AMD = 1 x 10~ 3 . The only periodic solutions of 
the secular Hamiltonian correspond to the oscillating Aw. 

Figure 7. Same as Fig. 3, but for AMD = 8 x 10 -3 . The domain of the secular resonance, where 
the angle Aw in libration around 0° and uj\ in circulation, are marked by 5. 
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Figure 8. Dynamical maps obtained for 7 mut = 5° (top) and J mut = 15° (bottom). The curves are the 
location of analytically calculated periodic orbits with Am = 0. The lighter regions indicate regular 
oscillation of Am (around or 180°), whereas the darker regions indicate its regular circulation. 
The domains, where planetary collisions occur within 0.5 Myr, are hatched. 

Figure 9. Same as in Fig. 8, but for I mn t = 45°. The blue curves are the location of analytically 
calculated periodic orbits with uj\ = 0, while the red curves of orbits with Am = 0. 

Figure 10. Dynamical maps of the region around the v Andromedae system on the (02, e2)-plane of 
initial semi-major axis and eccentricity of the outer planet d. Left: Map calculated with planetary 
masses of the current edge-on system, mi = 1.83Mj up and vrt2 = 3.97Mj up ; Right: the same, but 
with masses m\ = 2 x 1.83Mj up and m,2 = 2 x 3.97Mj up . The high-order mean-motion resonances 
are labeled on the top of the graph. The domains, where planetary collisions occur within 0.5 Myr, 
are hatched. The position of the edge-on outer v Andromedae system is indicated by a star. 



27 




28 




e^os 2to L 



4E-4 - 



3E-4 - 



OE+0 




~\ 1 1 

-0.5 -0.4 -0.3 -0.2 

e^os 20^ 




0.0 



Figure 2: Michtchenko et al. 
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Figure 3: Michtchenko et al. 
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Figure 4: Michtchenko et al. 
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Figure 5: Michtchenko et al. 
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Figure 7: Michtchenko et al. 
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Figure 8: Michtchenko et al. 
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